import os
import pandas as pd
from datetime import datetime
import numpy as np
import plotly.graph_objects as go
import plotly.express as px
import plotly.io as pio
from gpx_file_reader import GPXFile
# this allows rendering of graphics when exported to HTML
pio.renderers.default = "plotly_mimetype+notebook"
Get the GPX file we need, and convert it to a gpxFile object using the parser.
GPX_FILE = '20210808_181500_amtrak_dc_to_nyc.gpx'
# print("abspath: ", abspath)
# dirname = os.path.dirname(abspath)
# print("dirname: ", dirname)
# os.chdir(dirname)
# files = os.listdir()
# print("files:\n", files)
gpx_file_path = os.path.abspath('20210808_181500_amtrak_dc_to_nyc.gpx')
gpxFile = GPXFile(gpx_file_path)
gpxFile.print_info()
gpxDF = gpxFile.get_gpx_dataframe()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 8684 entries, 0 to 8683
Data columns (total 4 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 lat 8684 non-null float64
1 lon 8684 non-null float64
2 ele 8684 non-null float64
3 time 8684 non-null datetime64[ns]
dtypes: datetime64[ns](1), float64(3)
memory usage: 271.5 KB
preview:
lat lon ele time
0 38.962250 -76.859985 24.902832 2021-08-08 22:15:01
1 38.962430 -76.859764 25.251038 2021-08-08 22:15:02
2 38.962654 -76.859510 24.751282 2021-08-08 22:15:03
3 38.962833 -76.859310 26.062988 2021-08-08 22:15:04
4 38.963060 -76.859055 29.826477 2021-08-08 22:15:05
5 38.963220 -76.858850 26.231506 2021-08-08 22:15:06
6 38.963436 -76.858610 22.677307 2021-08-08 22:15:07
7 38.963654 -76.858380 22.543457 2021-08-08 22:15:08
8 38.963852 -76.858154 22.253967 2021-08-08 22:15:09
9 38.964058 -76.857910 20.727356 2021-08-08 22:15:10
Plot elevation vs lat,lon
fig = go.Figure(data=go.Scattergeo(
lon = gpxDF['lon'],
lat = gpxDF['lat'],
text = gpxDF['ele'],
mode = 'markers',
line = dict(width = 0.5,color = 'blue')
))
fig.update_layout(
title = 'Amtrak train ride',
geo_scope='usa',
width=1000,
height=800,
)
fig.show()
Proof that much of the elevation was negative. (Inaccurate data)
fig = px.line_3d(gpxDF, x="lon", y="lat", z="ele")
fig.add_trace(go.Scatter3d(
x = [gpxDF['lon'].iloc[0]],
y = [gpxDF['lat'].iloc[0]],
z = [gpxDF['ele'].iloc[0]],
name = 'DC'
))
fig.add_trace(go.Scatter3d(
x = [gpxDF['lon'].iloc[-1]],
y = [gpxDF['lat'].iloc[-1]],
z = [gpxDF['ele'].iloc[-1]],
name = 'NYC'
))
fig.update_layout(
title = "Elevation vs lat,lon",
width=1000,
height=800,
)
fig.show()
Plotting again, only taking into account positive elevation.
gpxDFElevationAbove0 = gpxDF[gpxDF['ele'] >= 0]
fig = px.line_3d(gpxDFElevationAbove0, x="lon", y="lat", z="ele")
fig.add_trace(go.Scatter3d(
x = [gpxDFElevationAbove0['lon'].iloc[0]],
y = [gpxDFElevationAbove0['lat'].iloc[0]],
z = [gpxDFElevationAbove0['ele'].iloc[0]],
name = 'DC'
))
fig.add_trace(go.Scatter3d(
x = [gpxDFElevationAbove0['lon'].iloc[-1]],
y = [gpxDFElevationAbove0['lat'].iloc[-1]],
z = [gpxDFElevationAbove0['ele'].iloc[-1]],
name = 'NYC'
))
fig.update_layout(
title = "Elevation vs lat,lon",
width=1000,
height=800,
)
fig.show()
# haversine formula: takes degrees
def dist_between(lat1, lon1, lat2, lon2):
R = 6371e3 #metres
phi1 = lat1 * np.pi/180
phi2 = lat2 * np.pi/180
deltaPhi = (lat2-lat1) * np.pi/180
deltaLambda = (lon2-lon1) * np.pi/180
a = np.sin(deltaPhi/2)**2 + np.cos(phi1) * np.cos(phi2) * (np.sin(deltaLambda/2)**2)
c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a))
return R * c # meters
def test_dist_between():
lat1 = gpxDF['lat'].iloc[0]
lon1 = gpxDF['lon'].iloc[0]
lat2 = gpxDF['lat'].iloc[-1]
lon2 = gpxDF['lon'].iloc[-1]
d = dist_between(lat1, lon1, lat2, lon2)
print("dist (meters): ", d)
dMiles = d/1609.34
print("dist (miles): ", dMiles)
#test_dist_between()
# velocity calculation
gpxDF["deltaDistMeters"] = dist_between(gpxDF['lat'].shift(), gpxDF['lon'].shift(),
gpxDF.loc[1:, 'lat'], gpxDF.loc[1:, 'lon'])
gpxDF["deltaTimeSeconds"] = (gpxDF.loc[1:, 'time'] - gpxDF['time'].shift()).apply(lambda row: row.total_seconds())
gpxDF["velocityMetersPerSecond"] = gpxDF["deltaDistMeters"] / gpxDF["deltaTimeSeconds"]
gpxDF["velocityMilesPerHour"] = gpxDF["velocityMetersPerSecond"] * (3600.0 / 1609.34)
gpxDFSpeedFiltered = gpxDF.loc[(gpxDF["velocityMilesPerHour"] >= 0) & (gpxDF["velocityMilesPerHour"] < 150)]
gpxDFSpeedFiltered.describe()
| lat | lon | ele | deltaDistMeters | deltaTimeSeconds | velocityMetersPerSecond | velocityMilesPerHour | |
|---|---|---|---|---|---|---|---|
| count | 8661.000000 | 8661.000000 | 8661.000000 | 8661.000000 | 8661.000000 | 8661.000000 | 8661.000000 |
| mean | 39.787519 | -75.563328 | -14.288275 | 35.543707 | 1.206905 | 34.525746 | 77.232086 |
| std | 0.432882 | 0.779921 | 10.459655 | 25.603769 | 5.413281 | 14.309065 | 32.008547 |
| min | 38.962430 | -76.859764 | -80.380710 | 0.222390 | 1.000000 | 0.131104 | 0.293272 |
| 25% | 39.401394 | -76.326996 | -22.057251 | 22.881986 | 1.000000 | 22.592893 | 50.538988 |
| 50% | 39.815760 | -75.432526 | -15.501739 | 37.802853 | 1.000000 | 37.569645 | 84.041111 |
| 75% | 40.094870 | -74.897545 | -6.697693 | 46.452316 | 1.000000 | 46.197891 | 103.341996 |
| max | 40.584797 | -74.303720 | 29.826477 | 1351.466293 | 303.000000 | 63.574951 | 142.213468 |
# Plot velocity
fig = px.line_3d(gpxDFSpeedFiltered, x="lon", y="lat", z="velocityMilesPerHour")
fig.add_trace(go.Scatter3d(
x = [gpxDFSpeedFiltered['lon'].iloc[0]],
y = [gpxDFSpeedFiltered['lat'].iloc[0]],
z = [gpxDFSpeedFiltered['velocityMilesPerHour'].iloc[0]],
name = 'DC'
))
fig.add_trace(go.Scatter3d(
x = [gpxDFSpeedFiltered['lon'].iloc[-1]],
y = [gpxDFSpeedFiltered['lat'].iloc[-1]],
z = [gpxDFSpeedFiltered['velocityMilesPerHour'].iloc[-1]],
name = 'NYC'
))
velocityColors = []
for i in range(0, len(gpxDFSpeedFiltered)):
v = gpxDFSpeedFiltered['velocityMilesPerHour'].iloc[i]
if v <= 50:
color = 'red'
elif 51 < v <= 100:
color = 'yellow'
elif 101 < v <= 150:
color = 'green'
else:
color = 'black'
velocityColors.append(color)
fig.add_trace(go.Scatter3d(
x = gpxDFSpeedFiltered['lon'],
y = gpxDFSpeedFiltered['lat'],
z = [0 for n in range(len(gpxDFSpeedFiltered))],
name = 'Route',
line=dict(color=velocityColors, width=0)
))
fig.update_layout(
title = "Velocity vs lat,lon",
width=1000,
height=800,
)
fig.show()
# Distribution of velocities with histogram
fig = px.histogram(gpxDFSpeedFiltered, x="velocityMilesPerHour", nbins=20)
fig.show()
https://github.com/NetTopologySuite/NetTopologySuite.IO.GPX/issues/34 This means that the GPX file I downloaded from the app actually does not conform fully to the schema.